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Abstract 

A detailed description of the numerical formulation of Version 2 of the ARIES/GEOS 
“dynamical core 1 ' is presented. This code is a nearly “plug-compatible dynamics for use 
in atmospheric general circulation models (GCMs). It is a finite difference model on a 
staggered latitude-longitude C-grid. It uses second-order differences for all terms except 
the advection of vorticity by the rotational part of the flow, which is done at fourth- 
order accuracy. This dynamical core is currently being used in the climate (ARIES) 
and data assimilation (GEOS) GCMs at Goddard. 


iii 



NTENTsONA 




eUCBm PAGE BLANK NOT FILMED 




Contents 


List of Figures 


vii 


1 Introduction 

2 Continuous Equations 

3 Conservation Properties 

3.1 Global Mass Integrals 

3.2 Total Energy 

3.3 The Vertically Integrated Pressure Gradient Force 

4 Horizontal Coordinate Systems 

4.1 Latitude-Longitude Coordinates 

4.2 Equations in Component Form 

4.3 Coordinate Rotations 

5 Vertical Differencing 

5.1 The Vertical Grid 

5.2 The Discrete Equations 

5.3 Energy Conservation and the Discrete u 

5.4 The Vertically Integrated Pressure Gradient Force 

6 Horizontal Differencing 

6.1 The Horizontal Grid 

6.2 The Continuity Equation 

6.3 The Thermodynamic Equation 


1 

3 

5 

5 

5 

7 

7 

7 

8 
9 

11 

11 

11 

15 

17 

18 

18 

21 

21 


v 


PRECEDING PAGE BLANK NOT FILMED 


PAGE ^ INTENTIONALLY BLANK 



6.4 The Tracer Equations 22 

6.5 The Momentum Equation 23 

6.5.1 The discrete equations in component form 23 

6.5.2 Fourth-order vorticity advection: the choice of a,/?, 25 

6.5.3 Energy conservation and the momentum equation 28 

6.5.4 Instabilities of the discretized vector- invariant form 32 

7 Polar Filters 35 

8 Time Differencing 36 

8.1 The Brown-Campana Scheme 37 

9 Using Subroutine DYCORE 38 

9.1 The Argument List 38 

9.2 Comments 39 

10 Final Remarks 40 


vi 



List of Figures 


1 The geographical and computational coordinate axes 10 

2 The vertical grid 12 

3 Position of the variables on the C-grid and the indexing convention 19 

4 Stencils of rjs used to define the Of, /?, etc 28 

5 Stencil showing the indexing of the of, /3, 7 , 6 factors, which appear in the 

energy conserving form of the inertial terms 31 

6 Stencil showing the indexing of the v and /x factors, which appear in the 

energy conserving form of the inertial terms 32 


vii 




1 Introduction 


The need to use modular coding techniques in atmospheric general circulation models 
(GCMs) has been recognized for some time. The two main benefits of this approach are 
that it. makes codes much easier to maintain and modify, particularly when the work must 
be done by several people, and that it makes it. easier to exchange and compare codes. Mod- 
ular coding practices are thus being used to some extent by most, atmospheric modelers; but 
most of the effort has gone into developing modular codes for the physical parameteriza- 
tions, such as radiation and cumulus convection. Kalnay et, al (1989) proposed a set of rules 
for coding physical parameterizations aimed at facilitating the exchange of codes. If these 
“plug-compatible" rules are followed in coding both the GCMs and the parameterizations, 
codes can he “unplugged” from one model and “plugged' into another with little effort. 

In this report we introduce the notion of a “dynamical core, which attempts to extend the 
ideas of plug-compatible parameterizations to the coding of the dynamics. Since the manner 
in which the equations of motion are discretized affects so fundamentally the organization 
of both the code and the data in a GCM, it is not practical to make the dynamics codes 
truly plug-compatible. They can, however, be made sufficiently modular that once the 
GCM’s data structure is modified to accommodate them the computational routines can 
be easily ported. 

The dynamical core we will be describing consists of a set, of subroutines that compute 
the time tendencies of winds, temperatures, surface pressure, and an arbitraiy number of 
tracers. The core is invoked by calling a single subroutine, and all data is passed to the 
core through that subroutine's argument list. Inputs to the core include: the dimensions of 
the grid, an array defining the vertical distribution of the layers, some physical parameters, 
the state variables at. two time levels, the time step, and the time tendencies of each state 
variable due to other processes. The core updates the time t.edencies to include the effects 
of the dynamics. It does not update the state variables, or perform any temporal or spatial 
damping that may be needed to control non-linear computational instability. 

Ideally one would want the core to be independent of the time differencing, but this is 
not possible except for the simplest explicit differencing. By passing the time step, two 
time levels of state variables, and tendencies due to other processes, the core can compute 
tendencies for other time differencing methods, such as the semi-implicit scheme commonly 
used in spectral models. The version we describe here is for explicit, time differencing only, 
but it can produce “economical explicit' tendencies (Brown and Campana 1978) for use 
with a leap-frog scheme. 

All state variables and tendencies in the argument list are on the core’s horizontal and 
vertical grids in our case, an unstaggered or Lorenz grid in standard sigma coordinates in 
the vertical and an Arakawa C-grid in the horizontal. Although the type of discretization 
is determined by the dynamical core and the GCM’s data structure must conform to it 
other things, such as the resolution or the placement of the layers in the vertical, need not 
be. In our dynamical core, these are controlled through the argument list and are thus 
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completely arbitrary (The resolution can even be altered from one call to another within 
a single program.) Thus, the core can be adapted to many applications by simply altering 
the call, rather than modifying the source and rebuilding the binaries. This should simplify 
the use of the core for those not interested in tinkering with the source code. On machines 
for which this dynamic allocation results in an inefficient code, this feature can be easily 
discarded and the array sizes fixed at compile time. 

The differencing used in the ARIES /GEOS core follows closely that which has been in use 
in the UCLA GCM since around 1980. Both the ARIES/GEOS core and the UCLA model 
are based on a C-grid, they use the same spherical grid, they both use Arakawa and Suarez 
(1983) vertical differencing, and they use similar forms of the momentum equation. There 
are, however, a number of important differences. The UCLA model, for example, uses a 
modified sigma coordinate, while the ARIES/GEOS core uses standard sigma. Rather than 
detailing these differences, we will present a complete description of the ARIES/GEOS core. 

This documentation is for Version 2 of the core. Version 1 was very similar, differing 
primarily in that it used a second-order scheme for advection in the momentum equation. 
Version 1 has been used for a number of applications at Goddard. In particular, it is the 
dynamics of the GEOS-1 GCM (Takacs et al. 1994), which was used for the re-analysis of 
the 1985-1990 period recently completed at the Data Assimilation Office (Schubert et al. 
1993). 

Two additional features are new in this version of the core. First, the code has been 
generalized so that the latitude-longitude coordinates on which the horizontal grid is based 
need not coincide with geographical latitude and longitude. By displacing the pole of 
the computational coordinate from that of the geographical coordinate, the more uniform 
grid near the computational equator can be brought to other regions of interest. This 
device has already been used successfully to model the polar stratosphere. Second, the map 
factors on the grid have been generalized to allow grids that are not uniform in latitude 
or longitude. Latitudinal stretching of the grid, with or without pole rotation, is already 
implemented; longitudinal stretching requires modifying the polar filters, a change that will 
not be available until the next version. The hope is that with these modifications we can 
easily convert the global models into regional, high-resolution models. 

The presentation follows the general form and notation of Arakawa and Lamb (1977), re- 
peating their discussion of the continuous equations in sections 2 and 3. Latitude-longitude 
coordinates and rotations between the computational and geographical coordinates are dis- 
cussed in section 4. The vertical differencing is presented in section 5 and the horizontal 
differencing, including a new fourth-order scheme for advection of vorticity is presented in 
section G. Section 7 describes the polar filters and section 8 the interaction of the dynamical 
core with the time differencing. Section 9 describes the FORTRAN interface. 
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2 Continuous Equations 


We use a o coordinate defined by 

V ~ PT 

(7 — 

7T 


(i) 


where p is the pressure, 7r = p$ — PT ? PS is the surface pressure, and pt is a constant 
prescribed pressure at the top of the model atmosphere. With px — 0 this coordinate 
reduces to the conventional a coordinate proposed by Phillips (1957). The code is written 
in this more general o coordinate for backward compatibility with earlier versions of our 
GCMs; in the current versions, however, we take px — 0. 


With this vertical coordinate, the continuity equation becomes 


dir 

dt 


-V ff • (ttv) - 


d(7r a) 
da 


( 2 ) 


where v is the horizontal velocity vector. Integrating (2) and applying the boundary con- 
ditions a = 0 at p = px an d p = ps, wc obtain the forms used in the model: 


#7T 

~di 



• (7TV) d(7 


(3) 


and 


(7T(T) = 




V a ■ (ttv) da. 


(4) 


The equation of state for an ideal gas is a = RT/p. where a is the specific volume, T is 
the temperature, and R is the gas constant. The following alternative forms will be used 
below: 


P „dP cj (dP\ cj (dP\ 
a _ R 0-_ Cp 8—- — [ da )^- a [ dn )^ 


(5) 


where 6 = T/P is the potential temperature, P = (p/po)*i K — R/ C p-, c p the specific heat 
at constant pressure, and p 0 is a reference pressure. In obtaining the forms in (5) we have 
used ^ = k £ and the relation 


(%) 

\da a V dit ) a 


(6) 


For the time being we are ignoring virtual effects. 


The hydrostatic equation is 


dp 


- -a, 
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(7) 


whore is the geopotcntial. Using (5) this can bo written: 


fd*\ 

\da) n ™ 


“7T 


c p 0 

a 




From (7) wo obtain 



which is the form used by Arakawa and Suarez (1983) for vortical discretization. 


(8) 


The thermodynamic equation is written in flux form to facilitate the derivation of a 9- 
conserving differencing scheme: 


0(tt6») d(nvO) ttQ 

—^7— = -Va • (7TV0) + 

dt da c p P 


( 9 ) 


where Q is the diabatic heating per unit mass. 


In addition to the equations of motion, the core computes tendencies for an arbitrary number 
of atmospheric constituents, such as water vapor and ozone. These are also written in flux 


form: 




dt 


da 


( 10 ) 


where ry^ is the specific mass of the kth constituent, and is its source per unit mass 
of air. 


The momentum equation is written in “vector- invariant” form, as in Sadourny (1975) 
and Arakawa and Lamb (1981), to facilitate the derivation of an energy- and enstrophy- 
conserving differencing scheme: 

^ = — (/ + ()kx v - a— ^ — V,t($ + K) — c.p0 S7 a P — — -r— , (11) 

dt da 7r da 

= -T)k x (irv) - - V „{$ + K) ~ c p Q V7r- f^T’ ( i2 ) 

where / is the Coriolis parameter, k is the unit vector in the vertical, 

C = (Vrr X V) • k 


is the vertical component of the vorticity along a surfaces, 


(/ + <) 


is an “external” potential vorticity, 

A' = l(vv) 

is the kinetic energy per unit mass, g is the acceleration of gravity, and T is the horizontal 
frictional stress vector. 
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The advoctive form of the momentum equation is obtained by using the identity 


in (11): 


- V(v • v) = (v ■ V)v — (V x v) x v 


dv 

~dt 


<9v 

= -/ k x v - (v • V„)v - a- - c p 0 

da 



a OT 
7 r da 


( 13 ) 


3 Conservation Properties 


3.1 Global Mass Integrals 


Integrating (2) over the entire atmosphere we obtain a statement of conservation of total 
atmospheric mass, 

jti nds= °’ 

where S represents the surface of the sphere. Similarly integrating (9) and (10), we obtain 
equations for the global mass- weighted integrals of (9 and 


d_ 

dt 



da ds — 



da ds 


and 


f f da ds — [ f (nS^) da ds. 

dt Js Jo Js Jo 


3.2 Total Energy 


To discuss conservation of total energy, we first write an equation for the temperature. 
Multiplying (9) by c p P we obtain, after some manipulation: 


d(nc p T) 


dt 


= -V a • (n\c p T) 
dP 


+ nc p 6 


dt 


„ „ . fdP 

+ v ■ V a P + a ( 


d(7vdc p T) 
da 

(w 


)J~ 


Q. 


(14) 


The first two terms on the right-hand side of (14) vanish when integrated over the entire 
domain. The first term on the second line represents the conversion to kinetic energy. Using 
(6). this term can be written as: 




Oft \ 


+ 7T <7 . 


(15) 



The thermodynamic equation can then be written in the more conventional form: 


<9(7 rCpT) _ 


-V ff ■ (irvcpT) - 


d(ir ac p T) 


+ (7 rua) + 7 tQ, 


( 16 ) 


dt ' '' ' 0C7 

where a; is the vertical p velocity: 

■^ + V * V 7T J , (17) 

and a is given by (5). 

Next we obtain the kinetic energy equation by dotting (7rv) with (12). The first term on 
the right-hand side of (12) is an “inertial force,” which does no work since it acts normal 
to the motion. The remaining terms give: 


d(n K) 
dt 


= tvK) 


d(n aK) 
da 


+ K 


( 


dn d(na)' 

_ + v„ • (»V) + -L-i 


- (.v) ■ (v„* + c„9 (|£) V.) - «v • 


(18) 


The last term on the first line of the right-hand side vanishes from conservation of mass 
(2). Terms on the second line represent the work done by the pressure-gradient force the 
conversion from potential energy. We next consider the first of these, — (7rv * V a <£): 


x / x f dn d(nd)\ 

~ (ttv • V a $) = -V<J * (7TV<I>) - $ H 

_ d(na$) dn .9$ 


d(n&<&) d(cr$) dn d$ 

da dt da 


da 


(™ + HI) 


= -V* ■ (t rv$) 

_ / * x t)(7r<7$) d(a$) dn 7 xc p 6 ( dP\ ( dir 


where we have used (2), (7), and the identity 

, 3(<t$) d<f> 

The kinetic energy equation (18) can now be rewritten as 

(7TV$) 

oa 

d(a$) dn 


(19) 


( 20 ) 


-V, • (t 


dt 


da 

cj /dp\ 


da dt 


n - 


a 


(^lh +ff (f +vV,r ). 


d(na<b) 

da 

+ V • V 7T 


f)V ■ 


dr 

da ‘ 


(21) 


Finally, using (17) and (5), we obtain a kinetic energy equation in which the conversion 
term is written in the same form as in (16): 


£>(7riQ 

dt 


= -V, 

d{a$>) dir 


d(ndK) 

( nvK ) V n 




(7TV$) 


c?(7ra<I>) 

da 


, , dT 

da dt ~ ™ ~ ~da' 


( 22 ) 


6 



Summing (16) and (22) and integrating over the entire domain, we obtain an equation for 
the rate of change of the total energy of the atmosphere: 


d_ 

dt 


[ [\c p T + K + $ s )—ds = 
Js Jo 9 




7r da 

ds y 

9 


where is the surface geopotential. 


(23) 


3.3 The Vertically Integrated Pressure Gradient Force 


The vertical integral of the pressure gradient force is irrotational in the absence of spatial 
variations of the surface geopotential. In pressure coordinates, 


[V S _ 

/ V p <f> dp = V 
JpT 


V 


= V 


[P s 

/ <3> dp 

JPT 

[PS 

/ ($ - $s) dp 

dp T 

[ n($-®s)d< 
Jo 


- <2>sVps, 

+ (p s ~Pt) V3>Si 


+ 7rV$S- 


(24) 


In the last two forms, the curl of the first term is zero and the second term vanishes if <J>s 
is constant. 

To satisfy this property in the discrete model, we will need to obtain this form starting from 
the pressure gradient force written in terms of gradients along a surfaces. 


J' 


7rV p <J> da — 


f 

V 


+ c p e (^-) VttI da, 


V 

V 


J (7T$) da — V7T J ^<I> + 

- <J> S V7T, 




da , 


f (7r<I>) 0 

Jo 

/ 7t(4>-^>s) 

L Jo 


da 


+ 7rV<3>s^ 


(25) 


where we have used (7) and (20). 


4 Horizontal Coordinate Systems 

4.1 Latitude-Longitude Coordinates 

In section 2, the continuous equations were written in vector notation. In this section we 
rewrite the equations as used in the code, in component notation in spherical coordinates. 
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Let (j) be the latitude and A the longitude. The differential elements of arc length along 
latitude circles and meridians are: 


dx = acos(f)d\ (26) 

dy = a d.(j), (27) 

where a is the radius of the globe. The differential element of surface area on the sphere is: 


ds = a 2 cos (j) dX d(fi. 


The gradient of a scalar field, t/h is 

A cty , $ 

= 7TT + “ ol’ 

a cos (p aX a 0(p 


(28) 


(29) 


where A and <f> are the horizontal unit vectors in longitude and latitude. 


4.2 Equations in Component Form 


Let u and v denote the zonal and meridional components of the horizontal velocity v in 
the (A, (j)) coordinate system. To write the scalar conservation equations for mass, 0 , and 
constituent concentrations in component form, we require only the form of the horizontal 
divergence of their mass-weighted fluxes: 


V 


ix\q 


1 / diruq c?(cos (j) irvq) 


a cos <j) \ dX 


+ 


d<j) 


)• 


(30) 


whore q can bo replaced by 1 for use in (2), by 8 for use in (9), and by q W for use in (10). 


The zonal and meridional components of the momentum equation (11) are: 


du 

dt 

dv 

~di. 


(/ + <)« 


du 


1 


o- 


\d($ + K) 


da a cos cjj 


dX 


.dv 1 \d{* + K) 


+ CpO 


* (f) £i 

\ dir )„ 3AJ 
(dP\ dn] _ 

V dn )„ d(j> J 


_ g_dT 
7r da 1 
q_dT 
7r da 1 


(31) 

(32) 


whei'e / = 2fisin^> is the Coriolis parameter and fi is the rotation rate of the coordinate 
system relative to the inertial frame. The relative vorticity, is 


. _ l f dv d(cos<fiu)\ 

( = V ff xv = ( — — ) . 

a cos 0 \aA c/0 / 


(33) 


The quantity (/ + C) regarded as the absolute vorticity, since / is the curl of solid-body 
rotation at rate ih 


f = V x (Qacos<f>X) 


d(Qacos 2 (j>) 


a cos <j) d (j) 


2fisin (j>. 


(34) 
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Eqs (31) and (32) can be written in the more conventional advective form by combining the 
relative vorticity and K terms: 


C v - 


1 OK 
a cos (j) dX 


-(u 


I dK 

a defj 


v ( dv 


a cos <j) \ dX 


<9(cos cf>u) \ 

04 > ) 


1 f du dv \ 

^V‘aA + ’’wJ 


v d(cos<f>u) u On 
a cos (/) d<f> a cos <f> d\ 

tan (f> u du v du 

a a cos (j) d\ a, d<f> 

) 

tan cj) 2 u & v v dv 
a a cos cf) dX a d(f) 


u 


a cos (j> 


dv 

d\ 


<9(cos (j>u) 

d4> 


\ If du 
) ~ a 


+ v 


dv 

d(f> 


4.3 Coordinate Rotations 

In the code, we allow the spherical coordinate system used for the computations to be dis- 
placed relative to the geographical latitude-longitude coordinates. Let (A, </>) denote longi- 
tude and latitude on the computational coordinate system and let (A, (j)) denote geographical 
longitude and latitude. The relation between the two coordinate systems is fully determined 
by specifying the coordinates of the geographical north pole in the computational system, 
which we denote by (Anp,</>np)i and by a third parameter, Aq, which represents a rotation 
about the geographical pole. The relation between geographical and computational axes is 
shown schematically in figure 1. Using these parameters, the two systems are related by: 


sin (j) = cos </>np cos </ J cos(A — Anp) + sin </>np sin (35a) 

cos(A + Ao)cos0 = sin</>NP cos</>cos(A — Anp) - cos </; NP sin (35b) 

sin(A T Ao) cos </> = cos</>sin(A — Anp), (35c) 

or 

sin (f) — cos ^np cos 0 cos(A T Aq — tt) + sin </>fvjp sin 0, (36a) 

cos(A + 7r - Anp) cost/; = sin^NP cos</;cos(A + Aq - tt) - cos</>np sin</>, (36b) 

sin(A + 7r — Anp) cos </> = cost/;sin(A + Ao — 7r). (36 c) 

Note that these transformations can bo inverted by making the exchanges: 

(A , </>) <=> (A,t/>), (37a) 

Anp <==> tt — Aq* (37b) 


The code assumes that all input variables are on the computational coordinate system, and 
that the input u and v wind components are also along the computational coordinates. It 
makes the same assumptions for the output variables. With these assumptions, the only 
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z 



Figure 1: The geographical and computational coordinate axes. 


effect on the code of displacing the geographical pole away from the computational pole is 
on the form of the Coriolis parameter: 

/ = 2f2sin</> = 2Q (cos</>np cos^6cos(A — A^p) + sin^^p sin (f>). (38) 

Although there is no other reference to quantities in geographical coordinates within the 
dynamical core, it is generally necessary to transform between the two grids when initializing 
a run or recording results in geographical coordinates. In addition to interpolation, this 
requires transforming the velocity components between the two systems. If we let x denote 
the local angle between the computational and geographical coordinates, the two sets of 
velocity components are related by 

u = cos x « — sin x v , 
v = sin x w + cos x w, 

and 

u — cos x ^ + sin x (40a) 


(39a) 

(39b) 
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cos xv — sin 


(40b) 


The angle x may be obtained from the relations: 

cos 0 cos x = cos <f) sin^>NP — c°s(A — Anp) cos */>np sin </*> (41a) 

cos^sinx — cos ^>np sin(A — Anp)- (41b) 


A subroutine package that implements these transformations is being developed and will 
be documented in a future report in the series. 


5 Vertical Differencing 

5.1 The Vertical Grid 

We use a Lorenz grid in the vertical, with both winds and temperatures defined at the 
same levels. The atmosphere between a = 0 and a = 1 is divided into LM layers by 
LM - 1 surfaces of specified constant <r t , l = 2,LM, so that a x = 0 corresponds to the top 
of the model atmosphere, and <7 lm + i = 1 corresponds to the lower boundary. At the LM 
layers we define the velocity, v/, the potential temperature, 0/, and the specific masses of 
all trace constituents, q\ k) . The vertical cr- velocity, o u is defined at the interfaces between 
the layers and at the top and bottom surface. This arrangement of the variables is shown 
schematically in figure 2. 

Note that both layers and interface levels are numbered from top to bottom and that, to 
maintain as close a correspondence with the code as possible, both are indexed by integers, 
rather than using the more common convention of assigning the interfaces half-integer val- 
ues. Level a/ thus corresponds to the top of layer l. Where there might be some confusion as 
to whether a quantity is defined at the layers or at the interfaces, we will mark the interface 
quantities with a hat (e.g., P). Vertical differences will thus be written as 


at the layers, and 

(8a)i = a/+i — a; 

at the interfaces. 

II 

£> 

1 

T 


5.2 The Discrete Equations 


Using (1), we define the pressures at the top, bottom, and interface levels as: 

pi = p T + 7 req, for l — 1, LM + 1. (42) 

1 The forms of the equations used in the code will be denoted by boldface numbers. 
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LM-2 


u, V, 0, q 


LM-1 


u, v, 0 , q 


LM 


u, v, 0, q 


LM-2 

LM-1 

LM 

LM+I 


Figure 2: The vertical grid 


The values of the Exner function at the interfaces are obtained directly from (42): 


Pl = 



for l = 1, LM + 1, 


(43) 


which can be written 

Pi = vr K (— Y = n K V t (43a) 

\PoJ 

for px — 0* The latter form is much more economical to use since the first factor on the 
right-hand side is independent of l and the second factor, “P/, depends only on a . The 
code takes advantage of such economies when Pt = 0 and uses the general form only when 
prjy ^o. We will thus write the — 0 forms when appropriate. 


Layer temperatures 7} are defined from the potential temperatures a s 7} = 0/P/, where P/, 
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the Exner function at the layers, is given by 


Pi = 


1 + K 


6(pP) 


8p 


, for l = 1, LM. 


Jl 


This is the form proposed by Phillips (1974). For px = 0 this becomes 


Pl = 


7T 


pS (i + «) 

where, again, Vi depends only on a. 


8 (cr K+1 ) 


8a 


= * K Vu 


(44) 


(44a) 


The vertical finite differencing of the thermodynamic and hydrostatic equations and of the 
pressure gradient force follows Arakawa and Suarez (1983). In their scheme, the potential 
temperature “interpolated” to the interfaces is 


n + 1 = 


(Pi +i - P±iWi + 1 + (Pi + 1 - Pi)0i 
Pi + 1 - Pi 


for / = 1, LM - 1. 


(45) 


As we show below, to conserve energy this form has to be used in the hydrostatic equation 
for the thickness between layers and in the vertical advection term in the thermodynamic 
equation. When pj = 0, (45) becomes: 


6 


t+1 


Pi+i — A+i 
pi + 1 - Pi 


0i + 1 + 


Pi + 1 - Pi 
Pi+i - Pi 


0i, 


so that the interpolation depends only on a, not pressure. 


(46) 


The hydrostatic equation is 

= $s + Cj>#lm(.Plm+1 - -Plm), (47) 

= $l+l + c-p0l+l(Pl+l ~ Pi), 

= $ /+1 +c p [(P /+1 -P /+1 )0 /+ , + (P, + i -Pi)0,}, for / = 1, LM — 1. (48) 

For pi = 0 we write this as: 

$i,m = + c p0+u^ K {P lm+\ - ’Plm), (47a) 

4*/ = + c r' nK [{Pi+\ ~ Pl+\)0i+\ + (Pi+i ~ Pi)0{\- (48a) 


The continuity equation is 


ebr . f 8(na)\ 


The equation for the tendency of 7r is obtained by summing (49) over all layers: 


LM 


= - £ v ' (7rVf ) 
at t.=\ 


(49) 


(50) 
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and the vertical a-velocity at the interfaces is given by 


d'K 1 

(t rcr);+i = ~cr l+ { — ^ V • (ttv*) (bo) £ , 

ot e - 1 


for l — 1, LM — 1. 


(51) 


The thermodynamic equation is 


dt 


— V * (ttv/0/) — 


<5(7Tfj0)\ 7 \Qi 

bo J ; + c p P, ’ 


for l = 1, LM, 


(52) 


where 0; is given by (45) and Pi by (44). We note that using (45) as the interpolation of 0 for 
vertical advection gives up conservation of 0 2 or some other function of 0. This is the price 
Arakawa-Suarez pays to use a local hydrostatic equation. The scheme of Arakawa (1972) 
satisfies all conservation properties of Arakawa-Suarez plus conservation of a function of 9, 
but it uses a non-local hydrostatic equation, in which one of the geopotentials depends on 
a weighted sum of all the temperatures in the column. 

For advection of the trace constituents, we are currently using the square-conserving form: 


dt 


= - V • (7TV/9 ; ( * :) ) 


8(ira 


8a 


Ji 


+ 7T«S 


■(*) 


for l = 1, LM, 


(53) 


where 



for 1 = 1, LM — 1. 


(54) 


The form of the momentum equation is: 


d\[ 

dt 


= -rp k x (ttv,) - [d,(v, - v,_i) + d /+1 (v, + i - v,)] 

2 (bo)i 

- V($ + K)l - C P 0,~V 7T - | (^) , 1 = 1, LM > 


(55) 


where 


dPi o l+] {P t+l - P t ) + atiPt - Pt) 


for l — 1, LM, 


dn 7r(8a)i 

as may be verified by differentiating (44). With p t = 0, (56) becomes: 

dPi _ 7T* a ;+ iQP f+ i -Vi) + (n{Vi-Vi) 

dlT 7T (^ o)l 

where, again, the second factor depends only on a. 


(56) 


(56a) 
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5.3 Energy Conservation and the Discrete u 


The derivation and justification of the vertical differencing scheme presented above is given 
in Arakawa and Suarez (1983). Here we will merely show that these differential-difference 
equations are energy conserving. We will also obtain a discrete form of u from the energy 
conversion term. 


We start, as in section 2, by writing the thermodynamic equation in terms of the tempera- 
ture at the layers, T/ = jP/ 0/. Multiplying (52) by c p P{ we have: 


d(ir<: p T,) 

Ot 


-V • (nviCpTi) - P, 


6(nac p 6)\ 


^ + v,.VP, 


+ ttQ,, (57a) 


= -V ■ (nvic.pTi) - 


f>(ir<7c p T) 
8(7 


+ 


(fi(T)l 


( 7 r ^ r )/+i(T/+i - prfi+i) - (nv)i{Ti - mj\ 


„ OP, 

\dir 1 


+ c p e,«— 

— + V/ • V7T 

L eft J 

i 


-f 7 tQi for l = 1, LM. 


(57b) 


Here the quantity T has been introduced to put the equation in flux form; this allows us 
to separate the energy conversion term. Comparing (57b) with (16), we see that the boxed 
terms of (57) correspond to ( 7 nua)i. 


To form the kinetic energy equation we take the dot product of (55) and 7 rv/: 
dnKi „ , x , / \ ] 

In + v ' ( ’ rv,> + (— ) ,j 

[(7r(7) /+ iV/ +1 • v, - (7rrr)/V; • v,_,] 


dt 


-V ■ (7 rviKi) + K t 

_ 1 1 
2 (6a), 

dP ( 8T \ 

- V • (7TV/$/) + $/V • (ttv l ) - C P^l' K ~Q^ V l • ^7r _ 9*1 • (^r) 
= -V • (71 -v,K,) - \yTT [( 7r< ^)/+l v /+l • v / - ■ v/— 1 ] 


V • (7rv/<J>/) — <$/ 


2 (6a), 
di r 


dt 


+ 


8(tt&)\] . dPi „ f 8T\ 


8(7 


Now introducing the edge-level geopotentials, <£, and using the identities: 

'6(a$)\ cr/ + i(4»/ + i - $/) + a,(®, - $;) 


$/ = 


8a 


(6a), 


(58) 


(59) 
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which is a finite differenced analog of (20), and 




for l = 1,LM, (61) 

where we have also used the continuity equation (49). 

Comparing (61) with (22), we see that the boxed terms in (61) correspond to the conversion 
term, — (nu)a)i. Energy conservation requires that the forms of ( 7 rua)i in (57) and (61) be 
identical. Comparing quantities in the first box in each equation and noting that d\ = 
M + 1 — 0, we see 1-hat we must require that 

c p (T)+j - P/0 /+ i) = — $/+i, for l = 1,LM — 1, (62a) 

Cpifi-PiBi) = for l = 2, LM. (62b) 

Similarly, considering the second boxes in each equation we must require that 

<7 f+ i($H-i - $/) + q/($/ - $/) = _ (9/7r ^ ) for 1 = LM. (63) 

(6a) t dir 

Using (56), we sec that this last is satisfied if the half thicknesses arc given by the following 
hydrostatic relations: 

($/+, - $/) = -c p 6i(Pi + 1 - Pi), for l = 1, LM, (64a) 

($/ — 6/) = —Cp8i(Pi — Pi), for/ = 2,LM. (64b) 

It is straightforward to verify that these forms of the half thicknesses are consistent with 
(48), the hydrostatic relations between layers. Furthermore, substituting (64a, b) in the 
2(LM — 1) relations (62a, b), the latter can be written: 

c p (Pi + 1 — Pt)8[+ 1 = — 4>/+i, for l = 1,LM - 1, 

which are identical to (48), and 

c p f,= ^(Pi + Pi-i)6i + \($i + $i-i) -$i, for l = 2, LM, (65) 

z z, 
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which define the LM - 1 interface temperatures we introduced in writing (57). 

Using the half thicknesses (64a, b) wc can rewrite (61) as: 

= -V • (irviKi) - ^ [(7T(7)/ + iV/ + i • V/ - (7rd),v, • v/_ j] 
at l [ball 


. { 6(nd$)\ ( 6{(7®)\ Ott /W\ 

V ■ (»v,*,) - - (— J ~oi - w ■ (fo), 


- c r t),* \ l* d M**' - f z 3\ + W (* + y, . yd , (66) 

p 1 7 1(60)1 on \dt )\ 

where the boxed term again corresponds to — (nuja)i. Using (56), this term can be written 


(nuja)i 


n dP, (»*)/+, (fl+i - Pi) + (ira),{Pi - Pi) /Sir 

■ ‘*’ r ar . 1+l( A +1 - a) + »,<«-«) + U +V ' v * 


c p 9i dPi . /Sir 

= 7 r-? — 7TfT + fT/ I — + V/ • V7T 

a i ok L V at 

Here tildes denote the weighted averages: 


(7T<7,<7/) = 


(7rg,<r) <+ i(^ + i - Pi) + (ir(7 , <t)/(F; - Pi) 
Pl+l-Pl 


which may be thought of as interpolations of these “edge quantities to the layers. The term 
in brackets on the right-hand side of (67) is a discrete form of io\. Clearly, the corresponding 
ai is: 

CpQi dPi 
on = • 

o\ on 

Note that the u>i obtained from (67) is defined at the layers, not at the interface levels, 
where <7/ is defined. 

Multiplying (57) and (61) by (60)1, summing over all layers, and integrating over the hori- 
zontal, we obtain the vertically discrete form of the total energy equation: 

, . LM f LM r fST\ 1 

/ Y'n(c p T l + K, + *s)(Scr) l ds= / V UQcSV, - l—JUSo^ds. 

ot Js •'S l—\ t \ / U 


5.4 The Vertically Integrated Pressure Gradient Force 

The discrete form of the vertically integrated pressure gradient force is: 

LM LM r dPi 1 

2>(V p $),(6a), = ^ir V$, + c p 0,— Vir (6<r), 
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Using (63), the identity (59) can be rewritten: 




Substituting this in (69) gives: 
LM 

5>(v p *m&7), 

/=i 


-m, 


. , ap, 


= V 


= V 


LM 

L/=l 
LM 

ll=] 


— $sV7T 

+ 7rV$Si 


which is in the same form as (25). 


6 Horizontal Differencing 


6.1 The Horizontal Grid 


The computational domain consists of a longitude sector extending from pole to pole, with 
cyclic boundary conditions applied in longitude. The size of the sector is determined hy 
specifying the AA-fold symmetry on the sphere. If M = 1, the domain is the full globe. 
Smaller sectors (A f >1) can be specified only if the coordinate pole and the geographical 
pole coincide. 

The grid divides this domain into JM latitude belts of width (A j = 1, JM, and IM 
zonal sectors of width (AA),, i = 1,IM. We define the grid (A,-, <f>j) as follows: 


o 

II 

1 

to| s* 

1 

II 

1 + (A^>)j, for j — 1, JM; 

<f>3M = 

^ | CN 
II 

(70) 

II 

1 

II 

+ 

+ (AA),-, for i = 1, IM; 

A IM+1 - 

II 

(71) 


Obviously, the (AA), and ( Acj>)j must be specified so that: 

™ 2tt jm 

= Jf 311(1 = 7r - 

J=1 j = 1 


The prognostic variables are located on an Arakawa C-grid. The pressure, temperature, and 
all tracers arc located at the points (A,, for all i and for j = 1, JM - 1, which excludes 


the poles. We will refer to these as the “p-points.” The “?/-points,” at which the zonal 
wind components are defined, are located between the “p-points” and on the same latitude 
circles, while “?;-points” are located between “p-points” and on the same meridians. The 
vorticity is defined at the “(-points,” on the same latitude circles as v and on the same 
meridians as u. In addition, the two pole points are considered both v- and (-points. This 
arrangement is shown schematically in figure 3, which also shows the relation of the grid 
intervals to the staggered variables and the indexing convention. Note that ( AA)y is defined 


I 

Pi, 2 


Mi- 1,2 

«i, 2, 

\ 

/ 

NCi — 1,2 

<72 

\ 

/ 

_A ™ / 

«»y,i 


\ 

/ 

\ 

/ 

\ 

/ 

Ci 

V\ 

<P0 

SOUTH POLE 


<t> J 


M 


J- 


(AA), 


r ~Pi,j u ij P* + 1 ,j~ 

(A <f>)j Ci y j v i + lJ’ 



u i,j-\ Pr + \ J—]- 

i : i 


INTERIOR 


/ 

\ 

/ u i— 1 ,JM — 2 


Pi 

1 

, JM—2 

1 

NORTH POLE 


Figure 3: Position of the variables on the C-grid and the indexing convention. 


at the u-points and (A <f>)j at the v-points. Variables defined at p-points are indexed the 
same as the u - point variables half a grid interval to the east , the ?/-point variables half a 
grid interval to the south , and the (-point variables to the southeast With this convention, 
we define the averaging operators as: 


At the p-points 


At the u-points 



2 Mi— 1 )j 

= 

+ V j + l)i 


(P’)ij = 

2 (Pi Pt+1 ) j 

(c j kj = 

2 (0 Cj + 1 )i 


and similarly for u and v at the (-points, and p and ( at the ?;-points. The difference 
operators, (<5px), (<5,p), etc., are defined analogously. 
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We will designate the grid spacings in arc length in longitude and latitude as Ax and 
Ay, respectively. These will be defined at the four points, p, u, v, C, and denoted by the 
appropriate superscript; thus, (A u x)ij is the longitudinal grid spacing at u-point (i, j). For 
the interior points we take: 


(A u x),j 

= a cos (f>j(A\)i, 

for j 

= 1, JM — 1, 

(72a) 

(A p x)ij 

= (A u x‘)jj. 

for j 

= 1, JM — 1, 

(72b) 

(A v x)ij 

= (A Px J )ij, 

for j 

= 2, JM — 1, 

(72c) 

(A 1 '?/), 

= a(A<f>)j, 

for j 

= 1, JM, 

(73a) 

(A p y) 3 

= (A V)j. 

for j 

= 2, JM — 2, 

(73b) 

(A u y)j 

= (A p y)„ 

for j 

II 

c-< 

£ 

i 

H -1 

(73c) 


The intervals (A 1 ^),-^, (A^rr^jM, (A p :r) tj o, (A p x)ijM, (A u y)o, (A u y)jM, (A p y)o, (A p y)jM, 
(A v y)u and (A”y)jM do not appear in the equations. The remaining intervals are special 
cases at the poles: 

(A”x) M = 0, (74a) 

(A”x)i i jM = 0, (74b) 

(A Py) x = a[(A^), + i(A0)2], (74c) 

(A p j/)jm-i = a . [(A^)jm + ^ i]- ( 74< *) 

The grid area elements are: 


(A l)ij = (A^ J ) tJ (A p y)j, for j — 2, JM — 2, (75a) 

(A Dij = (A u x)ij(A ,l y)j, for j = 2, JM — 2, (75b) 

(Al) u = (A v x)i j(A v y)j, for j = 1, JM, (75c) 


(A^)ij — (Ap )i,j, 

and at the poles 

(Ap);,i = (A^),-,i 
(Ap)i,JM-l = (A^)i,JM-l 

(A?), 

(A^)jm 


for j — 2, JM — 1, 

(75d) 

(A p x) M (A"y),, 

(75e) 

(A p x), i jM-i(A 1, t/)jM, 

(75f) 

1 r i ™ l 

2 IM ^ Ap ^J ’ 

(75g) 

1 r i im i 

2 fM^ (Ap), ’ JM_1 • 

(75h) 
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6.2 The Continuity Equation 


Wc use a simple eentered second-order form of the continuity equation: 


dnjj 

dt 



where u* and v* are: 


u*j t = ( A u y)j(n'u[)ij , for j = 1, JM - 1, 

v* j t = (A v x)ij(T j vi)ij, for j — 1, JM. 


(76) 

(77) 

(78) 


Summing (76) we have: 



which are the forms used in the code. The boxed term in (79) is the dynamical tendency 
of 7r computed by the core. Although there should be no other contributions to the 7 r 
tendency, the code adds the dynamical tendency to the n increment. Increments of the 
surface pressure (if any) passed to the core are not included in the that appears in 

(80); however, they can enter the calculation in the case of economical explicit differencing 
described below. 

We note that from (74a, b), ?;*, = ?;* JM = 0. Using these and cyclic conditions in longitude, 
the global sum of the divergence vanishes, and we have: 

|E(«EE^( a ?)v=°- <81) 

/ ! j 


6.3 The Thermodynamic Equation 

We are using the simple second-order, square-conserving (in the horizontal) form of the 
potential temperature equation: 
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where the 6 are given by (45). The boxed quantity in (82) is the dynamical tendency added 
by the core to the 7r 6 increment. Note that it is a mass-weighted tendency. Taking the 
mass- and area-weighted sum over the domain gives: 


d_ 

dt 


£(*/*)/ £ £(»*/), j (Kkj = £(*/*)< £ £ 


Sp P 1 , 




(83) 




Multiplying (82) by (A pC p P/),j and summing over the horizontal we have: 


EE 

* j 


d(nc p Ti) dPi di r 

— ^ 




* J 


+EE%« i,j,/ + (v*8 EE 


* jr 


* J L 


c p F/ 


6i(7T<J0) 
8 1 (7 


/J 


i^pkr (84) 


UJ 


In obtaining the first term on the second line, we have used the relations: 


Y. M = 0, a,*r = y] 6,-a 1 , (85) 

i i i i i 


and similar relations for j with the polar boundary conditions properly accounted for. The 
second term on the second line of (84) does not involve the horizontal differencing and 
so can be manipulated like the analogous term in (57); then using (62a) and taking the 
mass- weighted sum over all layers we have: 


d 


\ 


ot £(«/*)/ £ Yxwcjtiu (Ap)j v? 

i i j 

= £(«/*)i £ £ (a^q,) . . + £(M! £ £(*<**, Ai) itJ 

i i j ' 3 i i j 

7r^+i(^/+l “ p l) + nM p t ~ Pi) 




dPi\ dn ((u*6jP), 1 + {VI ijP), 3 ) 


/ dPi \ dir 

\ d,7r ) dt 


(86) 




6.4 The Tracer Equations 


We are currently using the same horizontal advection scheme for tracers as for potential 
temperature advection: 


d{^q\ k) )ij 


dt 




6i{u* qW )i + 8j(v*q( h '> J )i 


J tj 


f qW) 

8ia 




+(*$%’ < 87 > 


where the q ^ are given by (54). Summing over the domain we have 


d_ 

dt 


£(M £ £( 7r 9/ w ).j (Aj)ij = £(M ££ (^ w ) . . (Ap)ij. 


( 88 ) 
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6.5 The Momentum Equation 


The differencing of the momentum equation is patterned after that of the UCLA model. 
For the vertical advection term and the pressure gradient force we use a straightforward 
energy-conserving form, and for the inertial terms, a fourth-order version of the form of the 
Sadourny scheme used by Burridge and Haseler (1977). To derive this fourth-order scheme, 
we followed a procedure similar to that used by Takano and Wurtele (1982) to obtain a 
fourth-order version of the Arakawa and Lamb (1981) scheme. The resulting scheme is 
somewhat simpler than that of Takano and Wurtele, but retains most of the conservation 
advantages of the Arakawa-Lamb scheme. 


6.5.1 The discrete equations in component form 


We begin by writing the equations in a general energy-conserving form similar to that used 
by Arakawa and Lamb (1981): 


duiji 1 r ^ # 

Qt ~ (A^)7J l 1 i,j 

+ VijU i j_i — + + l — j -1 J + f ^J U i-Uj\ l 

~ Wh, 7m, \ -<«)], .. 


(89a) 


(A u x)ij 


Si{$i + K t ) + f>nr 


dir 


~r*r(ir) * f° r j — i, jm — i, 
X J {^)ij \ W J ,,j,l 


O y ,ji If, 

qI — (A v y)j J-l + 7i :j u ij + ^i-i j 

\,j ~ ~ ( Pij v ij + 1 + Wj-ltfj-l], 


(89b) 


(A v y)j 


dPi j 

SjiQt + Kri + Cpdt-^- 6p r 


I ~ 7=J\ (t~^) > for j — 2, JM — 1. 

Jij (* J kj V ' i,j,i 


Here the a, /?, 7 , 6 , e, <p, v, and )i are linear combinations of neighboring potential vorticities. 
The v and //. terms were added to the original Arakawa-Lamb formulation by Takano and 
Wurtele to expand the stencil enough to be able to require fourth-order accuracy. The 
potential vorticity at interior points is defined as 


(A^)ij(/tj + C i,j,l) c . 

==r , for j = 2, JM — 1, 

(*-A2 )ij 


(90) 


where the relative vorticity is given by 
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1 


(91) 


CiJ,l = Txi t- [(± v v)j(6iv)w - («i(A“xti,))ij] , for j = 2 , JM - 1 . 


At, the two poles, the potential vorticity is given by: 

VU 


(AlMh + Cu) 


^JM,/ = 


_ (A^)jm(/jm + Cjm,/) 


(90a) 


(90b) 


where 


Cl,/ = 
Cjm,/ = 


IM 


-<^IM 

1 1 IM 


(91a) 

(91b) 


In the code, polar values of the vorticities are replicated at all longitudes. 

The discrete form of the kinetic energy, K , is presented below. 

The discrete form of the Coriolis parameter is obtained by noting that / is the planetary 
vorticity; that is, the part of the absolute vorticity due to the rotation of the reference frame 
fixed to the planet. The velocity components of the rotating frame along the geographical 
coordinates are: 

Up = Via cos (f>, v p - 0, 

and from (39), the components along the computational coordinates are: 

Up = Via cos ^ cos x, v p = fia cos i^sinx- 

On the grid, we first define u p and v p at the p-points (and at the poles) and then interpolate 
to the u- and u-points: 


(u p )i j = fia cos 4> cos X , for j = 1, JM - 1, 


and 


(v p )ij = cos (j> sin x , for 3 = 2, JM - 1. 


(92a) 

(92b) 


Using the same discretization as for relative vorticity, we substitute these in (91) to obtain: 

Via 


fij — 


(A?)i 


C 


A v y S;(cos(/>sinx) - 6 J (A , ‘xcos^cosx ) 


, for j = 2, JM — 1, 


1 j 
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Qa 


(A ?)ij L 


A v i/^ t (cos<^Np sin(A — Anp)) 

— 8 j ^A w x sin <^np cos < f > — cos </>np cos(A — Anp) sin(/)| j 


J *j 


where we have used (41b). The form used in the code is 

-Qa f 

fi J = ^2 — | sin ^NP 6j( A u x cos <t>) 

— cos 0 np f A i; y 5| sin(A — Anp) + cos(A — Anp) 8j( A 1z x sin </>)! > , (93) 

L J 

which takes into account that, in latitude-longitude coordinates, (A — Anp) is independent 
of j and 4> is independent of i. The right-hand side of (93) is a discrete form of 


Q f . <9cos 2 <j> 

sin0 N p — — cos0 N p 


COS (j) 


d<f> 


<9sin(A — Anp) , , x . <9cos (/> sin </> 
— + cos(A — Anp) 


OX 


d(p 


(94) 


which reduces to (38). At the poles, we have: 


IM 


f] = (A^IM^^^ 1 ’ 1 (Sin0NP C ° S(f>] 

— cos (j ) np cos(A, — Anp) sin <p \ ) j , 


(93a) 


Qa 1 IM 

/jM = 7T2V — 7T7 J2 f(A“*) ti jM-l (sin^Np cos <£jM 
(A^)jm IM " 1 


-l 


-COS0NP coa(A, — Anp) sin^>jM_i)] , (93b) 

which arc obtained by substituting (92a) in (91a,b). 


6.5.2 Fourth-order vorticity advection: the choice of a,(3,. . . 

The procedure followed by Takano and Wurtele to obtain a fourth-order version of the 
Arakawa-Lamb scheme was to 1) form an equation for the second-order difference form of 
the vorticity, 2) assume that the flow is non-divergent and replace the advecting velocities 
by second-order differences of the streamfunction, and 3) equate terms in the resulting 
non-divergent vorticity equations with corresponding terms in the fourth-order Jacobian 
derived by Arakawa (1966). This procedure gives a scheme that is fourth-order only in the 
advection of a second- order vorticity by the non- divergent part of the flow. We will follow 
the same procedure, but without requiring that in the shallow-water equations the potential 
enstrophy be conserved for a general flow. Rather, we will require enstrophy conservation 
only for the case of non-divergent flow, as in the Sadourny scheme. We will refer to the 
resulting scheme as a fourth-order Sadourny. 
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As pointed out by Arakawa and Lamb (1981), the second-order Sadourny scheme corre- 
sponds to 


a li = Tlfoj+l + Vi,j + W+lJ+l), for j = 1, JM - 2, (95a) 

Pfj = i 2 (Vi,j + 1 + Vi,j + »7«-l j+i), for j = 1, JM - 2, (95b) 

jfj = n( 7 ?»,j+ 1 + Vi,j + Vi-ij), for J = 2, JM - 1, (95c) 

Sfj = + Vi,j + Vi+ i,j), for 3 = 2, JM - 1, (95d) 


where the j ranges take into account that v* j = = 0, and so some of the values 

neighboring the poles arc not used. For the second-order Sadourny scheme, 


and 


lj = <pfj = 0, 


u ?j = t l 'lj = 0- 


(95c) 


(95f ) 


Requiring enstrophy conservation only for non-divcrgent flow allows us to eliminate the 
e 2 j and terms. Dropping these terms is very attractive because we have found them 
quite troublesome. In three-dimensional calculations with the GCM, we found that they are 
responsible for much of the problem our implementation of the Arakawa- Lamb scheme had 
near the poles. These problems were eliminated when we used the second-order Sadourny, 
and so we used that scheme for Version 1 of the ARIES/GEOS dynamical core. We note 
that the v and // terms, required for fourth-order accuracy, arc similar in that they involve u 
in the u equation and v in the v equation. Indeed we found a similar, though lesser, problem 
when these terms were included. This problem, however, can be partially corrected, as we 
will discuss at the end of section G.5.4. 

Following the procedure outlined above, Takano and Wurtele found that the following re- 
lations must be satisfied to obtain fourth-order accuracy: 

<*i.j = Ci+i J + jjg [lGr/tj + i + 8(777+1^+1 + 7 )ij) - 3 (r/;_i j+i + r/,j+ 2) 

+ ( T 7i+lj-l + Vi±2j) - {Vi- l.j + Vi,j + 1 + Vi± lj+2 + 2,j+l )] 

- + ¥>)i+i.j> for j = 1, JM — 2, (96a) 

Pij — ~Cij + — [ 16 ^, j+i + 8(77,-! j+i + rjij) - 3 (t 7 j+ i j+i + T)ij + 2) 

+ ( 7 /i-1J-l + Vi- 2j) — (Vi+\,j + ViJ-l + Vi—LJ+2 + Vi-2J+i _ )] 

for j = 1, JM — 2, (9Gb) 

7 i.j = Cij + L [lGfy j + 8(77, _i,j + 7/ij+i) - 3 (r/ ?J _i + 77,+! j) 

+ ( 7 /»'-2,j + 1 + ^-l,i+2 ) - ( V»- 2J + *7»-lj-l + Vi,j+2 + Vi+\,j +\ )] 


26 



(96c) 


+ 2 ( c + f° r 3 = 2, JM - 1, 

Ki = -Ct+lj + L [16 ij,j + 8(?7;j + i + 7?, + ij) - 3 (t?,_i j + 

+ ( 7 ?j+ IJ+2 + 7?i+2,j + l ) - ( Vi+2J + 7 /i+l j-1 + Vi J+2 + Vi- 1 ,j+ 1 )] 

+ for j = 2, JM — 1 , (9Cd) 



- 24 (f-+ij 

Vi— 1 ,j ) i 

for j 

= 1, JM, 

(97a) 

t l ij 

= iitoj-i 

“ Vi,j+\ 

for j 

II 

to 

€-t 

1 

►— 1 

(97b) 


whore the quantities C, c, and (j) may be any combination of r/s. From (97a), /v 7 -j = z' 7 jyi = 0, 
and so terms involving u* j and ?q J M in (89a) disappear. 

Takano and Wurtele, following the same procedure as Arakawa and Lamb, use the choice 
of C, c, and (f) to obtain enstrophy conservation for a general flow. We follow Sadourny by 
setting 

ti-j = (fiij = 0 (98) 

and thus abandoning exact potential enstrophy conservation for divergent flows. This leaves 
only C to be determined. Note that the choice of C will only affect Of, /?, 7 , and 8\ the 
quantities v and /i» are the same for all fourth-order schemes. 

We will choose C\ j fairly arbitrarily, trying only to make the scheme as simple as possible. 
Figure 4 shows a stencil of the r/s neighboring u z j and the corresponding O', /?, 7 , and 8. This 
stencil covers all the r/s that appear in equations (9Ga-d). For comparison, we also show the 
stencils for the second-order Sadourny, as well as the Arakawa-Lamb and Takano- Wurtele 
schemes. The second-order Sadourny uses only the vorticities within the dashed rectangle in 
computing a:, /?, 7 , and 8 at (i, j). The Arakawa-Lamb and Takano- Wurtele schemes use the 
vorticities within the hexagonal region only two more points than the Sadourny scheme. 
Equations (96a-d) involve an additional eight points: the four enclosed in solid squares, 
which appear underlined in (9Ga,d), and the four enclosed in dashed squares, which appear 
underlined in (96b, c). 

We first require that C be such that the scheme uses only the Arakawa-Lamb stencil. This 
is satisfied if 

Cj j = CiJ + — [(?7i-l,j— 1 d" Vi — 2,j) — (^7* — 1 ,j+2 d~ Vi- 2 J+l) 

d" {ViJ+2 d~ Vi+ l,j+l) — d" ty'-hlj)] 1 (09) 

where Cij depends only on vorticities within the desired stencil. We then choose 

Qj — 24 bj Vi,j+ 1 ) “ (Vi- lj+i d" Vi,j)\ > (100) 
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Figure 4 : Stencils of 77 s used to define the ck, /?, etc. 


for which the scheme can be written in the particularly simple form: 
(a,p,y,f>)ij = |(a,/?,7»0?j ~ 5 («,/*. 7. 

whore 

nfj = —(*/»+! j + f/ij+2 + ty-l,j+i)> forj = l,JM — 2, 

/ ?£• = i(Wj+2 + »?i-ij+W+ij+i). for j = 1, JM — 2, 

7,j' = + *7f+l j + *7.-1 j+i). for J = 2, JM - 1, 

= ^(*7i+i j+i + *7.j-i + */.’— 1 j)i for 3 = 2, JM - 1. 


( 101 ) 

(102a) 

(102b) 

(102c) 

(102d) 


6.5.3 Energy conservation and the momentum equation 


To form a kinetic energy equation we multiply (89a) by ( A u x uj ), j and (89b) by (A y)j v t jj ■ 


-A4 

n -dT 
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+ U 7 ■ , \aiiV, 


[ a ij V i+\ j + 1 




— u 






(A 2 ), j 1 

(fa) 


^2 [(*<*)/“/(“/ ~ «/-i) + (7r<r) /+I M/(M/ + ] - «,)] 


— ?i 




Li*M 


W/ + ^/) + c„«/-3r ^ 


"£7f, 

(97T 


for j = 1, JM — 1, 




-.9 






, r 

- ,; >ji 


(Ag),-j 1 

(iS;rr)/ 2 


/ / »,i w t+l j ~ /'.-i jw.'-i jj / 

rri 




(ir&YjVtivi - w/ _ j ) + (7T^)/ + ,w/(w/ +1 - w/)] 

, for j = 2, JM — 1. 


dP, 3 

8j($i + Ki) + c pQi-fa fa 




(103) 


(104) 


The global integral of the kinetic energy tendency at each level can now be obtained by 
adding (103) and (104) and summing over all % and j. Doing this, we have on the left-hand 
side: 


JM— 1 

: e 

t j = i 


al 2\ JM-1 / c) 12' 

E E +EE [ A ’*^r 


dt 

JM-1 




1 J 


A 5^“? + Aji»? 


■ ? £ -4 

Z J = 1 

J M — 1 f) TSC 

= EE 

« j=i 


1 j 


(105) 


since (A 2 ),,i = (A^) t> jM = 0. Here K c is the discrete form of the kinetic energy per unit 

1 


mass: 


Klhl (A J)y 


A?, ^«/ 2 + A 2 E 2 j , for j = 1, JM - 1. 


( 106 ) 


We can perform a similar manipulation under the summation of the right-hand side, since 
in the u equation all reordering is done over i : which is cyclic, and in the v equation, terms 


29 



at the poles that appear when the sums are reordered as averages at the p-points contain a 
factor of either v* or which are both zero at j = 1 and j = JM. 


EE 


[I ~ ^ --EE» 

1 1 


« j 


A2 


u t 


\Sia), 


+ A ?, 


hw : 

- EiEj \ + (7T-T) ( + |A'^U,(m 1+] - Jl() J 

+(< 7T(7 ) l Alv l (v l — vt-iY + {na)i+iAlvi(vi + i - vrf) J 

+ E , Ej (*/ + Ki) («i«i + «i«D - (“? + ?M J )] . ’ ( 10? ) 


1 J 





+ Ei Ej 
+ Ei Ej 
-EiEj 

— E» Ej 


( Ap Kf)ij 


d*i,i . 1 / c * t c j* \ . f&li™)} 

■sr + (a^'*' + 6i "' + (— J w 


(#/ - Kf)ij(8iu1 4 - f>jvf)ij 
1 1 


(Sig)i 2 [(( 7r< ^* +1 A« u t u <+ 1 * ~ (ira)/A5 «/«/_, ) 

+ ^(vr(j)/ + iA2 u/Uj + i J - (7rd);A2 ) t; / u < _ 1 J ^j 


Ap^/ 




1 J 


(108) 


The second line of (108) drops out using the continuity equation, and the third line vanishes 
if we take 

K = K C . 

Note that the a, (3, 7, 8, v, and /x terms do not appear in (108). The a, (3, 7, 8 terms cancel 
between the u and v equations when summed over the globe, assuming an appropriate 
treatment at the poles. This cancellation can be clarified by considering the stencil shown 
in figures 5. 

In figure 5, the ot, f3, 7, and 8 are shown between the us and us they affect. Each term appears 
in the equation for the neighboring u multiplied by the neighboring v and in the neighboring 
v equation multiplied by the neighboring u, but with opposite sign. For example, (3 ,j 
appears in the u z ] equation multiplied by u*j + , , and in the u,j + i equation (obtained by 
incrementing j by 1 in (89b)) multiplied by u*j. (Note that the convention we use here 
and in the code is that the a, (3, 7, 6 are indexed the same as the u equation in which 
they appear.) The v and /x terms vanish when (103) and (104) are summed over i and j, 
respectively. The stencil for the v and // terms is shown in figure 6. Note that v and /x are 
defined at the ^-points. 
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Figure 5: Stencil showing the indexing of the a, /?, 7 , 8 factors, which appear 
in the energy conserving form of the inertial terms. 


Multiplying (108) by (fycr)/, summing over all layers, and using (76), we obtain the following 
form of the total kinetic energy equation: 


l * 3 

£<«££{( a’)*, [§ + (W + ?V) 


where V is the totoil kinetic energy dissipation: 


15 = E^'EE 

/ i j 



, (109) 


The first term in the airly brackets in (109) involves no horizontal differencing and can be 
manipulated in the same way as in the horizontally continuous case discussed in section 5.3, 
equations (61) through ( 66 ). Doing this gives 


^ E^)' E E( A ?)*J 7r <j (*/ + - E( 5 ^)' E E(* ^ A J)« j 

" l i j I ' i 

(7T<T)/ + l(Pf + l - P/) + (7rcr)f(P; - P<) 1 / dPA 

7 r( 6 i<r)i j V d 7 r / 



d7T U ; * (^,71 ’ + Vj Sj'K 3 

dt + 7T 


(HO) 


Notice that the factors in curly brackets in the conversion terms in (110) and ( 86 ) are not 
identical, and thus the horizontally discrete form of the equations does, not conserve energy 
exactly. The difference stems from the 6 V 7 T term in the pressure-gradient force. If 

this term were discretized as 

VtiP, 


rather than 
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Figure 6: Stencil showing the indexing of the v and fi factors, which appear 
in the energy conserving form of the inertial terms. 


conservation would be exact. We have not done this because the conserving form does not 
allows us to maintain the irrotationality constraint on the vertically integrated pressure- 
gradient force. In any case, we will be abandoning exact energy conservation in applying the 
polar filter and other modifications described in the following sections. The non-conserving 
form of the equations also allows us to simplify the implementation of the economical explicit 
scheme described below. 

Comparing (110) with (66), we see that the horizontally discrete form of (n ua)ijj analogous 
to (67) is 




( dPi 

\ d/ dir 


7TG + <7/ 


i-3 



V.J T Vj Sj 7T J 


7 r A 2 

" p 


ij 


(in) 


where (ira, are given by (68) evaluated at each grid point. This is the form used to 

evaluate the uj returned by the code. 


6.5.4 Instabilities of the discretized vector-invariant form 


Hoi lings wo rth- Kallberg instability. 

The discrete form of the equations presented above suffer from the computational instability 
discussed by Hollingsworth et al. (1983). This instability is associated with the horizontal 
differencing of the inertial terms and is common in schemes written in vector-invariant form. 
It manifests itself as noise with short meridional scales and long zonal scales in regions of 
strong zonal flow. 
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Hollingsworth ot al. (1983) present, a heuristic argument suggesting that a spurious source 
of momentum is producing the instability. To illustrate this argument we consider the one- 
dimensional shallow-water equations written in vector-invariant, form. Ignoring sphericity 
we have: 


!-</ + 0» = o, 

Ov Oh, 0 1 2 1 2 

m +if+0u = ~%~oi, ( 2 u + r >• 

Expanding terms in the v equation, we obtain the advective form: 


Ov Ou Oh Oku 2 



The argument of Hollingsworth et al. is that the underlined terms do not cancel in the 
discrete form of the equations, resulting in a spurious source of momentum. The following 
linear system illustrates the effect of having such a term in the equations: 


Ov. 

m~ fv 

Ov 

m +{H 


Oh 

Ot 


0 , 


-<] 


-H 


Oh 

Oy 

Ov 


u 


Oy 


On 
Oy 1 


whom the linearized error arising from non-cancellation is assumed to be of the form 
The dispersion relation for this system is: 


a 2 = f + t)Hl 2 + tfUl, (112) 

where a is the frequency, l is the wavenumber in ?/, and i = Note that unstable 

solutions occur when / and U are non-zero. 


The simplest way of obtaining cancellation is to modify the form of the kinetic energy that 
appears on the right-hand side of the momentum equation. Hollingsworth et al. present 
such a modification for the second-order Sadourny scheme. We will follow a similar approach 
that is currently being used in the UCLA model to remove the instability from the Arakawa- 
Lamb scheme. Rather than requiring exact cancellation, this approach requires only that 
the terms cancel when the equations are linearized in cartesian coordinates. Although a 
much weaker constraint than exact cancellation, this removes the instability with a smaller 
modification of the kinetic energy term, and so maintains better conservation. 


For one-dimensional flow, 


„ 1 Ou 2 

U( '~ 2 Oy ' 


We wish the discrete form of the inertial terms to satisfy this relation when linearized about 
a uniform zonal flow. In this case, the discrete form of u( at latitude j is: 


U (aq-i + fij - 1 + jj + 6j) = U + 0+1 + 0-i)> 


(113) 
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since 


«»-i = Pj - 1 


7 3 = 8 3 


^( 2 Cj + Cj-i) - ^(O-i +Cj + Cj+i), 

1 _5_ _J_ 

12^ _1 + 24^ J 24^ +1 ’ 

^(2Cj 4 0+0 - 24 ( 0—1 + C j + Cj+i). 

1 5 1 

12 Ci+I + 24° ~ 24 Cj '- 1 ' 


It is worth noting that the quantity (atj-\ + fij. 1 + 7 ? + O') is independent of the choice of 
C; thus that choice could not be used to optimize the form of K. Substituting 


0 — ( u j u j— 1 ) 

in (113), we have: 

U (rtj-i 4 /?j_i 4 7 j + O') = ~ u - uj. 1 ) 4 u j+ 1 - Uj 4- Uj 1 - Uj— 2 ) , 

= -[/ ^((u J+ i + lOuj + «j-i) - («j 4 10uj_i 4 Uj. 2 )). 

Since we wish to require that 

U (c*j - 1 4 ft-i 4 7j 4 O') = -(*> - ifj-i), 


we see that the linearized form of the kinetic energy must be 


Kj = U 


5 1 !/ , x 

+ g 2 U - ?-1 + U]+x ’ 


This is the linearization of 


where 




5 0 1.2 

<H + 6 U ? 


Uj = ~(uj_i 4u J+ i). 


Inserting the map factors and applying this scheme in both directions, we have: 




(114) 


where K r is given by (106) and 


Kjj = 


(A l) 


p/h3 


* 5 *? 


+ a; A.-.? 


for j ~ 2, JM — 2 , 


(115) 




where 


1 


These corrections are applied at the interior points, j = 2,JM — 2. At the poles the 
scheme for the inertial terms differs significantly and so a different modification of K is 
required. However, we have not derived such a scheme and, for simplicity, have been using 
the conserving form: 


“ Kl\J an< ^ KiJ M-l,/ = ( 118 ) 

It seems likely that at least some of the problems remaining at the poles are due to this 
crude treatment. 


The v term, at the poles . 


As discussed above, we found that the v term given by (97a) tends to produce zonal noise 
near the poles. This is similar to the problem we had encountered earlier with the e and 
4> terms in Arakawa-Lamb scheme, and seems to be related to the fact that these inertial 
terms, although conservative, are not acting normal to the flow. The problem, however, is 
confined to the poles, and the culprit seems to be the v factor at j — 2 and JM — 1. We 
arc circumventing this problem by replacing these us by interpolations between interior and 
pole values (note that v vanishes at the pole): 


^*,2 = J't, 3 


(A'y)i 


J'i.JM— 1 




(AP?y)l + {APy) 2 ' 

(A p j/)jm-i 


(A p y)jM-i + (A p t/)jM— 2 ' 


(119a) 

(119b) 


since Uj\ = 0. 


7 Polar Filters 


In the dynamical core, we apply a polar Fourier filter to the tendencies of all prognostic 
variables. The purpose of the polar filter is to avoid linear computational instability due to 
the convergence of the meridians near the poles. The filter acts poleward of 45° latitude, and 
its strength is gradually increased towards the pole by increasing the number of affected 
zonal wavenumbers and the amount by which they are damped. The polar filter is also 
applied to the diagnostic cj. 

Let iftij denote a single level of any of the tendencies to be filtered. Its zonal Fourier 
expansion is: 

IM/2 

ipi,j = ^2 Vvnj cxp(-imA t ), 
m=0 

where is the complex amplitude of the mth zonal wavenumber and % = y/—l. The 

filtered tendency, rfij, is 

IM/2 

rfij = y GXp(”T77fcAj), (120) 

m=0 
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where the filter coefficients are given by: 


( co s(4>j) i y 

y cos (<^>450) sin(mj^j-) J 


( 121 ) 


The rationale for the form of the filtering coefficients is the same as in Arakawa and Lamb 
(1977). Crudely speaking, for n = 1 this filtering “slows-down” the zonal propagation of 
each zonal component enough to satisfy the same stability criterion as is required by the 
shortest wave at a particular latitude in this case 45°. Arakawa and Lamb (1977) use the 
meridional grid size, rather than the zonal grid size at 45° latitude, and they use the filter 
more selectively, treating only the gravity waves by smoothing only the mass divergence and 
the pressure-gradient term, rather than the entire tendencies. We found that applying it to 
the entire tendencies made the model more stable without adversely affecting the results. 
Also, if the filter is applied more selectively, one must be careful of how it interacts with the 
inertial instabilities discussed in section G.5.4. We have also found it useful to take n = 2, 
which provides additional filtering near the pole and seems to enhance the stability. 

Note that in this form the polar filter is only appropriate for grids with equal spacing 
in longitude. Modifying the filter for more general grids should be straightforward; for 
example, one can first interpolate in the zonal direction to a uniform grid as fine as the 
finest spacing on the non-uniform grid and then apply the same scheme. Another possibility 
is to design a local filter that approximates the Fourier filter for a uniform grid. We are 
currently developing such schemes. 


8 Time Differencing 


As discussed in the introduction, the task of the dynamical core is limited to computing 
the dynamical contributions to the tendencies of all state variables. With time continuous, 
these depend only on the instantaneous state. For some explicit time differencing schemes, 
such as the leap-frog br the Euler-backward, these tendencies depend on the state at a 
single time level and are independent of the particular time-differencing chosen. But for 
multi-level explicit schemes and for all implicit schemes, the tendencies depend on the time 
differencing and on the size of the time step. 

In the ARIES-GEOS core, the argument list includes two time levels of the state variables 
and the time step. This argument list can accommodate the most commonly used, explicit 
or implicit, time-differencing scheme. In Version 2 we support only explicit schemes. When 
the specified time step is less than or equal to zero, the core returns the explicit, single-level 
tendencies. When a positive time step is specified, the core returns “economical explicit” 
tendencies, using the scheme of Brown and Campana (1978), as described below. We plan 
to include semi-implicit differencing in future versions of the core. 
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8.1 The Brown-Carapana Scheme 


Brown and Campana (1978) proposed an explicit, scheme to use in conjunction with leap- 
frog differencing which relaxes somewhat the instability condition for gravity waves. The 
scheme does not affect the advection or inertial terms. The idea is to average the pressure 
gradient force over three time levels. This can be accomplished very simply when using a 
standard leap-frog scheme by updating the mass field (potential temperatures and surface 
pressure in our case) first, then computing the pressure gradient at the three levels and 
taking an appropriate average. An even simpler strategy, which is equivalent, when the 
pressure gradient force is linearized, is to average the three time levels of the mass field and 
only compute the pressure gradient once, using the averaged values. We take the latter 
approach, simplifying it still further by using the averaged mass field only for selected parts 
of the pressure gradient, calculation. 

We finite difference the momentum equation in time as follows: 


v " + 1 V ' = ] V 7 r* + ( all other terms )", 

At V <T/r 


( 122 ) 


where the <[>“ are obtained from the hydrostatic equation: 






<f> s + c ; Am*(J\m+i - 

$,”+1 + c p [{Pi + 1 - Pi + \) n 0i+]' + {Pi + 1 - PlY'Ol ] 5 


(123) 

for l- 1,LM - 1, (124) 


and 


= n l . 


(125) 


Here the ( / denotes the three time level average 

7 f‘ = a E (n n+i + K n ~ 1 ) + (1 - 2a F ,)ir n . (126) 

As pointed out by Brown and Campana (1978), the averaging parameter, a E , must, be 
chosen carefully to maximize the time step. The best value depends on the strength of the 
Robert (1966) time filter used to control the splitting of solutions with the leap-frog. We use 
a E = .2475, which is near optimal for a time-filtering parameter of .05. We emphasize that 
the Robert, filter is not done by the dynamical core, only the Brown-Campana smoothing 
of the pressure gradient that enters the momentum tendencies returned by the core. 

Another point to note is that the dynamical core will use a completely updated mass field 
at time level n + 1 only if the tendencies due to all other (non-dynamical) processes are 
input. However, it is not essential to do this to obtain the time step improvement, since 
what matters is the treatment of the gravity waves. 
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9 Using Subroutine DYCORE 


The dynamical core is accessed by a call to FORTRAN subroutine DYCORE. Inputs to 
the routine include the state variables, various parameters that specify the horizontal and 
vertical grids, and some physical constants. The outputs of the routine are the updated 
time tendencies of all state variables. The state variables themselves are not modified. 

A typical GCM application would use the following procedure at each time step: 

• Compute the “physics’ - tendencies, storing them in UOI, VOI, POI, QOI (PH would 
normally be zero, since sources or sinks of mass are usually ignored). 

• Call DYCORE, which would add the dynamical contributions to the tendencies. 

• Add the contribution of any horizontal diffusion. 

• Perform the time stepping, including any time filtering that may be needed, to obtain 
updated values of the state variables. 


9.1 The Argument List 

The following is a listing of the SUBROUTINE declaration and the “banner” that is included 
at the top of the DYCORE source code. The banner describes the argument list and the 
routine’s storage requirements. 


SUBROUTINE DYCORE ( 


IM, JM 

, JDIM 

, LM ,SIG jPTOP j 

P KM,DT 

OMEGA 

, CP, 

RGAS , 

AE, 


PHS 

, PKH , 




PIB , 

,U0B 

, VOB 

, POB 

.QOB, 

PIM , 

,U0M 

, VOM 

, POM 

, QOM, 

PII, 

UOI , 

VOI, 

POI, 

QOI, 

OMG , 

VOR, 

DIAG 

) 



C INPUT ARGUMENT DESCRIPTION 


C 

IM ... 

. Number of Grid Intervals in Zonal 

Direction 

C 

JM ... 

Number of Grid Intervals in Meridional 

Direction 

C 

JDIM . 

. Meridional (Second) Dimension of Input 

Fields 

c 

LM . . . 

. Number of Vertical Levels 


c 

SIG . . 

(LM+1) : Sigma at Interfaces. SIG(1)=0; 

SIG (LM+1) =1 

c 

PTOP . 

Model Top Pressure at Sigma = 0.0 


c 

KM . . . 

. Number of Scalars, Including H20 , but not Theta. 

c 

DT ... 

Time-Step fron n-1 to n+1 (Seconds) 


c 

OMEGA. 

Rotation rate (rad/sec) 


c 

CP 

Specific heat at constant pressure (J/(kg K)) 

c 

RGAS. . 

Gas contant ( J/(kg K)) 


c 

AE 

’Earth’ radius (meters) 


c 

PHS . . . 

(IM.JDIM): Surface Geopotential (m * m/sec**2) 
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C PKH . . . (IM,JDIM,LM+1): (P/POO) **KAPPA (Not used if PTQP=0) 


c 

PIB 

(IM, JDIM) : 

Mass (Psurf-Ptop) 

mb 

at 

Current 

Time-Level 

c 

UOB . . . 

(IM, JDIM.LM) 

: Zonal Wind 

m/s 

at 

Current 

Time-Level 

c 

VOB 

(IM, JDIM.LM) 

: Meridional Wind 

m/ s 

at 

Current 

Time-Level 

c 

POB . . . 

(IM, JDIM.LM) 

: Potential Temp . 

K 

at 

Current 

Time-Level 

c 

QOB 

(IM, JDIM.LM, 

KM): Scalar Fields 


at 

Current 

Time-Level 

c 

PIM . . . 

(IM, JDIM) : 

Mass (Psurf-Ptop) 


at 

Previous 

Time-Level 

c 

UOM 

(IM, JDIM.LM) 

: Zonal Wind 


at 

Previous 

Time-Level 

c 

VOM 

(IM, JDIM.LM) 

: Meridional Wind 


at 

Previous 

Time-Level 

c 

POM 

(IM, JDIM.LM) 

: Potential Temperature 

at 

Previous 

Time-Level 

c 

QOM . . . 

(IM, JDIM.LM,: 

KM) : Scalar Fields 


at 

Previous 

Time-Level 


C OUTPUT ARGUMENT DESCRIPTION-- Tendencies are in per second. 

C PII ... (IM,JDIM): Updated Surface Pressure Time-Tendency 

C UOI ... (IM , JDIM , LM) : Updated Zonal Wind Time-Tendency 

C VOI . . . (IM , JDIM , LM) : Updated Meridional Wind Time-Tendency 

C POI ... (IM , JDIM ,LM) : Updated Pi-Weighted Theta Time-Tendency 

C QOI ... (IM , JDIM , LM , KM) : Updated Pi-Weighted Scalar Time-Tendency 

C OMG ... (IM , JDIM ,LM) : Omega Diagnostic (mb/sec) 

C VDR ... ( IM , JDIM ,LM) : Vorticity Diagnostic (1/sec) 

C DIAG . . Logical (On/Off) Flag for Diagnostics 


C NOTES: 

C (1) JDIM is to be used for DIMENSION Purposes ONLY 
C Vectorization will be performed over IM*JM. 

C (2) The Vertical Layers are numbered from T0P(1) to BOTTOM(LM) . 

C (3) All Time-Tendencies are INCREMENTED (bumped) . 

C The Momentum Time-Tendencies ARE NOT mass-weighted. 

C l*he Potential Temperature and Scalar Time-Tendencies ARE 

C mass-weighted (by PI) . 

C (4) JM is 180 degrees divided by the meridional grid size. 

C (5) UXX(I.J) are located half a grid interval EAST of PXX(I t J). 

C VXX(I f J) are located half a grid interval SOUTH of PXX(I t J). 

C (6) If PTOP>0, the PKH MUST be defined. 

C (7) The previous time level fields (PIM , UOM , etc) are used for the 
C economical explicit calculation done in conjunction with 

C leap-frog steps. If you are not doing leap-frog or do not 

C wish to have economical explicit tendencies, pass DT<=0. 

C SPACE REQUIREMENTS: 

C (1) Takes IM* JM*19+4*JM+2*IM words from the heap for STATIC storage; 
C these are kept throughout the run. 

C (2) Takes IM* JM* (LM+25) + 3*LM + 1 words from the heap for DYNAMIC 
C storage when PT0P=O; for PT0P!=0, add IM*JM*LM words. 

C All of this storage is freed before returning. 



9.2 Comments 

• The argument. JDIM is provided in case the calling program has the state variables 
organized by level rather than variable, or has extra latitudes. 
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• With PTOP sot, to /,oro, tho code economizes by raising only the surface pressure to 
the k. If PTOP is non- zero the code expects the (^) K at the interfaces to be passed 
in the array PKH. This is done to avoid costly recomputations of this factor, which is 
usually needed outside the dynamical core. If PTOP is zero, the input variable PKH 
is completely ignored; it need not be defined, or even allocated as an array. In this 
case we are not concerned about duplicating the exponentiation. 

• The code sets the grid according to the values of IM and JM. This is done on the first 
call or whenever IM or JM change from the previous call, and the grid is saved in a 
dynamically allocated static space. The vertical grid is inexpensive to compute and 
so it is reset on each call. 

• In the current version the time step is used only for computing Brown-Campana 
tendencies. It is also used as a flag for these tendencies, which are only computed 
when DT is greater than zero. If DT is less than or equal to zeio, the code returns an 
explicit tendency computed with the PIB, UOB, etc. state variables, and the PIM, 
UOM, etc. are ignored. 

• Access to the grid rotation or irregular grid features is not provided through the 
argument list, in this version. The location of the pole and the zonal symmetry can be 
changed by altering the appropriate parameters in subroutine GRIDH. An irregular 
grid in latitude can be specified by modifying the definition of the Tj . also in routine 
GRIDH. Irregular grids in longitude are not supported in this version. 

• A message passing version is also available. 


10 Final Remarks 


We have presented a detailed description of the Eulerian grid-point dynamics currently being 
used as the “dynamical core” of global atmospheric modeling at Goddard. This dynamical 
core is written in a modular form that allows easy implementation, assuming the GCM’s 
variables are already stored on an Arakawa C-grid. In addition to a detailed discussion of 
the numerics, including the introduction of a new fourth-order scheme, we have described 
the usage of the routine. 


The next version of the core will include: 


• A fourth-order scheme for the horizontal advoction of potential temperature and trac- 
ers 


• full telescoped grid in latitude and longitude 

• An improved treatment of the poles in the momentum equation 
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All of those changes are in fairly advanced stages of development. 

The FORTRAN source code is available over the network. To obtain access to it, or for 
further information, send e-mail to: 


M ax . J . Su arez(b)gsfc . n as a. gov or L awrence . L . TakacsOgsfc . n as a. gov 
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